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Introduction 

In this paper, we give an algorithm to generate connection formulas of 
generalized hypergeometric functions pFp_i for degenerated values of param- 
eters. We also show that these connection formulas give a fast method for 
numerical evaluation of generalized hypergeometric functions near oo. 

Several methods to evaluate generalized hypergeometric functions are 
known; see, e.g., the famous text book "Numerical Recipes" |5]. As Van Der 
Hoeven proved, evaluating series gives a fast method when the precision is 
big. Since we need high precision values of generalized hypergeometric func- 
tions, we will use series expansions of generalized hypergeometric functions 
near oo for numerical evaluation. We call the series expansion of pFp_i(z) at 
z = oo a. connection formula. Several methods to obtain connection formulas 
in degenerate cases are known among experts, but an algorithmic method 
which is fast and relevant for numerical evaluation is not known. Van Der 
Hoeven gave a method to construct series solutions around regular singular 
points by introducing an order among logarithmic monomials x™(logx)" [S]. 
Note that Saito, Sturmfels and Takayama found an analogous method which 
is generalized to several variable case [H Chapter 2] . By utilizing these meth- 
ods, we will give a new method to obtain connection formulas in degenerate 
cases. 

Methods discussed in this paper are used in our numerical checker for 
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a digital formula book project [0] and are implemented in by Risa/Asir. 
All timing data are taken by Risa/Asir on a machine with the following 
specification; AMD Athlon MP 1800+ 1533.40-MHz, memory 2Gb PC/AT 
machine. 



1 Hypergeometric Function 

The function defined by the following series is called the Gauss hypergeo- 
metric function. 



F(«,/3,7;^) = E 



-z 



1^0 (7)fc(l)fc 
{a)k = {a){a + l)---{a + k-l) 

The analytic continuation of this function is also called the Gauss hyperge- 
ometric function. This function satisfies the following differential equation. 

{6,{6, + 7 - 1) - z{6, + + (3, r,z) = 

oz 

This differential equation is called the Gauss hypergeometric differential 
equation. By expanding the products 6^ in terms of z and we obtain 

{z{z - 1)— + (7 - (« + /? + l)z)— - aP}F{a, P, z) = 

Hence, the Gauss hypergeometric function has singularities at 2; = 0,1, cxd. 
It is known that the hypergeometric function have the following integral 
representation (Euler integral representation) 

F{a, f3, 7; z) = ^ /' - tr—\l - tz)-^dt 

1 (aji (7 — a) Jo 

Our goal is numerical evaluation of the generalized hypergeometric func- 
tion defined by the following series 

,F,_,{au ■ ■ ■ , a„ A, ■ ■ ■ , /3; z) = f ^"^^^ " " " ^"^^^ 



,t^o(A)^---(/?.-i)fc(l) 
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The famous book "Numerical recipes" [5] says that "a fast, general routine for 
the complex hypergeometric function 2-P'i(a, b, c; z), is difficult or impossible". 
One difficulty is that the generalized hypergeometric series converges only 
in \z\ < 1, then the method of evaluating series can be used only when 
\z\ < 1. However, in case of the generalized hypergeometric functions, there 
are connection formulas, by which we can express the hypergeometric series 
pFp_i in terms of a set of series which converges at 2; = oo. For example, in 
case of 2-^1, the connection formula is as follows. 

Proposition 1 (connection formula, see, e.g., [T] ) Assume a — (3, (3 — 
a,7 ^ Z<o. Then, we have 

F{a,(3,-f;z) = 

+ 

The condition a — (3, (3 — a, j ^ Z<o means that a — (3, (3 — a, j are not in 
the set {0, —1, —2, —3, . . .}. When the condition is satisfied, we say that the 
parameters are generic and when the condition is not satisfied, we say that 
the parameters are non-generic or degenerated. 

We can use the connection formula to evaluate numerically the hyper- 
geometric function near z = cxd in generic case. Degenerated cases will be 
discussed in the rest of this paper. We note that the evaluating hypergeo- 
metric series can be accelerated by the binary splitting algorithm [2], [3]. 

2 Connection Formulas of pFp_i 

In this section, we study connection formulas of pFp_i between and 00 and 
numerical evaluation by using the formula. In the case that parameters are 
generic, these formulas are well-known. They can be obtained by using the 
Barns integral representation (see, e.g., pj Chapter2 4.6]). When parame- 
ters are non-generic, there is no complete list of connection formulas nor an 
algorithm to obtain connection formulas in the literatures. We will give an 
algorithm to derive connection formulas when parameters are non-generic. 
The Algorithms 1 and 2 seem to be implicitly used among experts to study 
global behaviors of generalized hypergeometric functions, but Algorithm 3 
will be new and it gives a fast routine. 



I(2)I|z|f(a.a-.+ l.a-3 + l)(-.)-" 

r 7-/5 r a 
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Let us review the connection formula in the generic case. It will be the 
starting data to generate connection formulas in non-generic case. 

Proposition 2 (connection formulas of pFp_i, see, e.g., pLj) Assume ai — 
o^jl^ 7^ j) Z<o. Then, we have 
pFp^i{oii, ■ ■ ■ , Oip, Pi, - ■ ■ , Pp-i, z) 

ai-ai + l,---,ai- ai-i + - a^+i + 1, ■ ■ ■ , - + 1; 
Algorithm 1 The case of ai — 0:2 G Z: 

1. Use contiguity relations to make ai = 02- 

2. Multiply (0:2 — «i) to the both sides of the connection formula in the 
generic case. 

3. Replace (0:2 — «i)r(— ai+a2) by V{—ai+a2 + l) and (02 — cn)r(a;i— 0:2) 
by -r(ai - ^2 + !)• 

4. Apply ^ for the both sides and take the limit 02 — > ai- 
This method can be generalized to more degenerated case. 

Algorithm 2 We assume ai = a2 = ■ ■ ■ = oiq {q < p) and ^i^ di — aj^Z 

iq<i,J <p,i7^ j): 

1. Multiply Iljllll'j^i_^_i{aj — a^) to the both sides of the connection for- 
mula in the generic case. 

2. Replace (oj— aj)r(— aj+ctj) by T{—ai+aj + l) and {aj—ai)T{—aj+ai) 
by —T{—aj + + 1). 

3. Apply f for the both sides and take the limit 02 — * ai, ■ ■ ■ , «g — 

i—l i 

Note that the left hand side of the output of the algorithm is 

(nfrizl) pFp_i(ai, . . . , ai, Oq+i, . . . , Op, . . . , /3p_i; z) 

and hence the right hand side of the output gives a series expansion of this 
function around z = 00. 
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We conjecture that a repetition of applying the Algorithm [2] and of ap- 
plying contiguity relations yields a connection formula for any degenerated 
case. 

This method requires a complicated symbolic differentiation, but com- 
puter algebra systems are good at it. One can say that our method of using 
connection formula of pFp_i to evaluate numerical values around z = oo is a 
hybrid computation of symbolic computation and numerical computation. 

Example 1 The following is a connection formula in a degenerated case 
which is obtained by the Algorithmic] 

2Fi{ai, ai, Pi, z) 

lim - — 2-ri(ai, 1 + «! - Pi, 1 + "1 - ^2, 



r(Q;i)r(— «1 + OL2^ai Qci^' 

lim - — 2-Fi(a2, 1 + "2 - A, 1 - «! + "2, 



r(ai)r(— ai + (3i) "2^01 da2 

7r(A) + r(A)^(ai) 



r(ai)r(-ai+/5i) 

7r(/5i) - r(/3i) iog(-^) + r(/?i)7/;(-«i + i3i 



^)-"^2i^i(«i,l + «i-A,l,lA) 

(-z)-"^2Fi(«i,l + «i-/3i,l,lA). 



r(Q;i)r(-ai + i3i 

Here, ^^{z) is the derivative of log(r(z)) and 7 is — ^/'(l). 

Our computer experiments show that Algorithm [2] is not efficient. The 
next algorithm is an efficient version of deriving connection formulas in the 
degenerate case. Algorithm [2] is used only to get c°'s in the following algo- 
rithm. 

Algorithm 3 In the case of ai = 02 = ■ ■ ■ = aq{whereq = p), series solu- 
tions at z = 00 can be written as 

00 

= (-;.)-- ^{cj, + cl log(-z) + ■ ■ ■ + c;_i \og^-\-z)}z-^ 

i=0 

1. Derive recurrence relations for d- with respect to i by a generalized 

hypergeometric differential equation (see, e.g., J^, J^). c^- (i > 0) are 

,0 



determined by 's. 



2. Obtain by applying a part of Algorith 



m\ 



3. Get other coefficients c* by using the recurrence relations for c* with 
respect to i. 
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3 Examples of Using Connection Formulas in 
Degenerate Cases 



Mathematica implementation for numerical evaluations of generalized hyper- 
geometric functions is known to be very nice. Here we compare numerical 
evaluation by our formula (derived by our Algorithm 3) , evaluation by Math- 
ematica and evaluation by numerical integration. 

Example 2 We try to evaluate the value of 2^1(10/3, 10/3, 7/2, 13+137^). 
We will present timing data of evaluation by the series expansion derived by 
Algorithm 3 and Mathematica. We present timing data on Mathematica only 
for reader's convinience; it is nonsense to compare our timing data by our 
algorithms and those by Mathematica, because an algorithm used by Math- 
ematica is not known and implementations are done on different languages. 



terms 


time 


value 


5 


0.0019SOC 


().()0()0464654506S423618485+0.000098886376836B4298440V-1 


10 


0.002762sec 


0.00004646537447334307802263261624 
+0.00009888640350652418669794640828^^ 


20 


0.00417Bsec 


0.000046465374473393490391242220236585714989 

+0.000098886403506421825123991664023061 171848 \/-l 


40 


0.008503sec 


0.00004646537447339349039124213865727073014,'')885()337(i()3(i(i()133374 
+0.0000988864035064218251239916232587199904578128942387359317473v/^ 


80 


0.02377sec 


0.0000464653744733934903912421386572707301458850337603824784541 

+0. 0000988864035064218251239916232587199904578128942387442282741 v^-l 



Here, terms mean the truncation degree of hypergeometric scries. The 
precision is set to terms + 10. This timing data includes time to set c* in 
Algorithm 3. As we see in the table, the computation is done in less than 
0.03 second. 

Timing data by Mathematica 4.0: 



pn TL^ioii 


thiic 


xaluc 


IIJ 


I). 1)4 Sc^eoiul 


U.l)l)l)l)4l)4().';4 — ().()0()0<ISSS(i4v' -i 


25 


0.04 Second 


().()()()()4646537447339349039124214 
+0. 00009888640350642182512399162 


50 


0.07 Second 


0.000046465374473393490391242138657270730145885033760382 

+ 0.000098886403506421825123991623258719990457812894238744V-1 



We also try to evaluate the same value by using Barns integral represen- 
tation and adaptive rule. 
Timing Data: 



precition 


time 


value 


4 


0.3885SCC 


0.00004509188619331184102+0.00009303289997220799852 


5 


0.466sec 


0.00004508711196227677087+0.00009302906708707449476^^ 


6 


0.702sec 


0.00004508788940142053260+0.00009302899442986751422 
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Figure 1: a timing comparison of different algoritlims 
"datalr" : Algorithm 3 
"datalm" : Mathematica 

"datalb" : Euler's integral representation and trapezoidal rule 



Example 3 Wetry to evaluate the valueof2-Fi(7/2, 7/2, 31/5, 1.3+1.8v^), 
which will be more difficult than the previous example for series expansion, 
since 1.3 + 1.8^/^ is closer to the boundary of the domain of convergence. 



terms 


time 


value 


5 


0.001905sec 


-0. 3879786816479458591-0. 27(i75435;-i84fl()170:-i(i8v'-l 


10 


0.00272SOC 


-0.3770255218705445491-0.2823972087891714305%/ -1 


20 


0.004093SOC 


-0. 3709544095052939938707251207-0. 2822863971357611098403957229v'-l 


40 


0.008422sec 


-0.37695442761307946514230306490910664462 
-0.28228642179392542114229797464872838012 V-l 


80 


0.023Blsec 


-0.376954427613081226577499361640669979083664967652834206588 
-0.282286421793927502415734929810558926710399341428368162640^-1 



The precision is set to terms +10. The convergence is slower than the 
case of z = 13 + 13^—1- We also try to evaluate the same value by using 
Euler's integral representation and trapezoidal rule. 
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Figure 2: a timing comparison of different algoritlims 
"data2r" : Algorithm 3 
"data2m" : Mathematica 

"data2b" : Euler's integral representation and trapezoidal rule 



sample size 


time 


value 


2000 


13. Usee 


-0.3769544276724114820-0. 2822864217813210745 


4000 


85.41sec 


-0.3769544276222648211-0. 2822864217919839264 V~ 


8000 


346.4sec 


-0.3769544276144992000-0.28228642179362812147" 



The precision is set to 19. The numerical integration requires more CPU 
time, but the accuracy seems to be better than the series expansion. This 
data tell us that series expansion at 2; = oo should not be used around \z\ — 1. 

Timing data by Mathematica 4.0: 



precision 


time 


value 


10 


0.01 Second 


-0.376954 - 0.282286V-1 


25 


0.01 Second 


-0.3769544276130812265774994 -0.2822864217939275024157349V - 1 


60 


0.01 Second 


-0.37695442761308122657749936166305176029442543627034 
-0.28228642179392750241573492983143989797240125238853 V-1 



8 



The Figure below is a visualization of the difference |r20 — rlO|. Here, r20 
is the truncation of our connection formula of 2-Fi(7/2, 7/2, 31/5, x + y\/^) 
at the degree 20 and rlO is that at the degree 10. 




The difference is close to zero in \z\ » 1, but it is larger in the neighborhood 
of 1 2; I = 1. Developing a method for high-precision evaluation around = 1 
will be a future problem. 

Example 4 We will evaluate values of 3^^2(7/2, 7/2, 7/2, 31/5, 36/7; z). We 
compare the adaptive integration method of the Barns integral and numerical 
evaluation of our degenerated connection formulas derived by Algorithm [31 
Numerical integration (precision is set to 19): 





time 


value 


z = 130 + 130V- 1 


0.2219sec 


0.000001881452796232078012 
0.000006893851655427774880y~ 


Z = 13 + l3^/^ 


0.1605sec 


0.007312359138710618527 
-0.006400306129230969169^/" 


z = 1.3 + 1.3^/^ 


0.112sec 


-1.097506492885595820 
+0.6369234717153367928^/^ 



Series expansion (terms are set to 20, and precision is set to 19): 
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time 


value 


z = 130 + 130v^ 


O.OlSrisec 


0.00001345106300346753915 
0.000006796099418228839164v^ 


z = 13 + 13v^ 


0.01404sec 


0.007350815068974895610 
-0.006282360701607166085v^ 


z = 1.3 + 1.3v^ 


0.004741sec 


-1.097992622097576377 
+0.6364759787999937697v^ 


Timing data by Mathematica 4.0 (precision is set to 10): 






time 


value 




z = 130 + 130v^ 


0.1 Second 


0.0000134511 
+6.79615 lO-^v^ 




z = 13 + 13v^ 


0.07 Second 


0.00735089 

- 0.00628229v^ 




;z = 1.3 + 1.3v^ 


0.07 Second 


-1.09725 
+0.636973^/^ 





We conclude that our formulas derived by algorithm 3 give a fast method 
to evaluate generalized hypergeo metric functions around oo. 
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